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SECTION  I 
INTRODUCTION 

A new  finite  element  supersonic  kernel  function  method  in  lifting  surface 
theory  was  presented  in  Ref.  1.  This  manual  contains  the  Finite  Element 
Supersonic  Kernel  Analysis  Program  (FESKAP),  developed  for  the  new  method. 
Descriptions  of  the  main  program  are  presented  as  well  as  on  the  preparation 
of  input  necessary  to  execute  the  program.  A sample  run  is  included  to  illus- 
trate the  usage  of  the  program.  Descriptions  of  each  subroutine  are  pre- 
sented in  Appendix  A,  and  the  program  listing  is  contained  in  Appendix  B. 

The  purpose  of  the  computer  program  is  to  generate  the  generalized 
force  coefficients  at  one  specified  Mach  number  and  reduced  frequency  for  a 
given  planform  and  a given  set  of  modal  deflections.  The  program  is  appli- 
cable to  any  isolated  arbitrary  planform  in  supersonic  flow  with  subsonic/ 
supersonic  leading/trailing  edges.  No  thickness  effect  is  accounted  for.  The 
unsteady  motion  is  assumed  to  be  harmonic  for  the  aralysis. 


1.  Young,  V.  Y.  C.,  and  M.  R.  Brashears , "A  New  Finite  Element  Supersonic 
Kernel  Function  Method  in  Lifting  Surface  Theory,"  LMSC-HREC  TR 
D496650,  Lockheed  Missiles  & Space  Company,  Huntsville,  Ala.,  December 
1975. 
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SECTION  II 

PROBLEM  DESCRIPTION 


According  to  Ref.  1,  the  finite  element  formulation  of  the  integral  equa- 
tion in  the  lifting  surface  theory  is  given  as 


8 


(e) 


( 3Tx  +ik)  £i(x-y)  = 47  Z*. 


, J-  v ct(e)  \!e} 


(i) 


where  k is  the  reduced  frequency,  L (x,  y)  is  the  i * modal  function  and 


is  the  column  vector  containing  the  nodal  lift  values  of  an  element,  due  to  a 

(e) 


.th 


unit  displacement  in  the  i"“  mode.  ^ denotes  summation  over  the  elements 
within  the  forward  Mach  cone. 


The  row  vector  containing  the  integrated  kernel  coefficients  is  defined 


as 


Ct(e)  = jj  N£  (x,y)  . K (x -xo,y -yo)  dA 


(2) 


(e) 


where  N^(x,  y)  is  the  row  vector  of  shape  functions  and  K is  the  kernal  function. 


The  generalized  force  coefficients  are 

N 


% ■ S //  f.  (x,  y)  Nt(x,  y)  dx  dy  . xj 


T(e) 


(3) 


i=1  A<e> 


N 

where  ^ denotes  summation  over  all  elements. 


i = l 


Equations  (1),  (2)  and  (3)  form  the  framework  for  the  computer  program 
development . 


SECTION  III 

PROGRAM  DESCRIPTION 


The  program  as  presently  set  up  is  extremely  compact.  For  example, 
for  a case  of  239  nodes  constituting  222  elements,  the  program  size  is  slightly 
under  20K  (decimal)  words.  Variable  dimensions  are  used  in  all  the  sub- 
routines, so  that  the  user  needs  only  to  change  the  first  dimension  statement 

in  the  main  program  to  fit  in  a new  planform.  No  overlay  nor  auxiliary  file 
is  used  in  this  program. 


For  convenience,  the  mode  shapes  are  built  into  the  program 
are  represented  by  the  set  of  1,  x,  x2,  y2,  x2y2,  y,  Xy  or  the  set  of  1,  x,  x:  , 

Y , y,  xy.  The  sets  are  identified  by  the  number  of  modes  they  contain.  The 


These 
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J Ain 

user  has  the  option  of  specifying  the  first  set  (NMODE  = 7)  or  the  second  set 
(NMODE  = 6) 


For  a given  planform  and  Mach  number,  the  user  must  first  define  the 
characteristic  mesh  that  best  fits  the  planform.  The  nodes  are  then  numbered 
starting  from  left  to  right  with  the  foremost  points  and  proceeding  downstream. 
The  elements  are  numbered  in  a similar  manner.  Fill-in  triangular  elements 
with  horizontal  sides  at  the  supersonic  trailing  edge  are  numbered  last  since 
they  cannot  be  an  influencing  element  to  any  collocation  point. 


Influencing  elements  within  the  forward  Mach  cone  of  a collocation  point 
are  determined  automatically  by  the  program.  The  element  containing  the 
collocation  point  as  its  most  downstream  node  is  defined  as  the  pivotal  ele- 
ment. All  elements  with  element  number  less  than  the  pivotal  element  number 
are  scanned.  Thus  all  candidates  are  either  forward  of  or  on  the  same  level  as 
the  pivotal  element.  Each  in  turn  is  further  tested  by  a logic  statement  to 
see  if  it  is  within  the  Mach  cone. 
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As  explained  in  Ref.  1,  a table  of  weighted  kernel  function  coefficients 
is  first  tabulated  for  later  table  look-up  during  the  solution  process.  The 
si/.e  of  this  table  is  governed  by  two  parameters  IMAX  and  JMAX,  IMAX  is 
the  maximum  number  of  characteristic  elements  in  the  chordwise  direction 
as  determined  from  the  mesh.  For  most  planforms,  this  is  given  by  the 
number  of  elements  on  the  centerline.  To  find  JMAX,  locate  the  most  ex- 
treme collocation  point,  which  is  usually  the  one  close  t,o  the  tip  of  the  trail- 
ing edge.  JMAX  is  the  number  of  layers  of  characteristic  elements  necessary 
to  cover  all  the  elements  within  the  forward  Mach  cone. 

To  take  advantage  of  the  symmetry  and  anti  - symmetry,  the  lift  calcu- 
lation is  performed  only  on  the  right  half  of  the  planform.  Each  node  has  its 
associated  mirror  image  with  respect  to  the  centerline.  Lift  values,  with 
positive  or  negative  sign  depending  on  whether  the  mode  is  symmetric  or 
anti-symmetric,  are  simply  substituted  in  for  the  mirror  node.  For  con- 
sistency, a node  on  the  centerline  has  itself  as  the  mirror  node. 

A schematic  flow  chart  of  the  program  is  shown  in  Fig.  1. 

Description  of  Variables 


BETA 

p =Vm2 - 1 

COEF 

Array  of  integrated  kernel  function  coefficients 

IANGLE 

Number  of  angle  of  element;  IANGLE  = 3 for 
triangle  and  IANGLE  =4  for  quadrilateral 

IBUF 

Buffer  array  for  printing  out  the  list  of  influ- 
encing element  number 

ICHECK 

Option  parameter  for  quick  mesh  check  run 

IMAX 

Maximum  number  of  regular  characteristic 
elements  in  the  chordwise  direction 

INFO 

Array  of  element  nodal  information 

JMAX 

Maximum  number  of  regular  characteristic 
elements  in  the  Mach  line  direction 

KMAX 

Number  of  element  information  data  cards  to 
be  read  in 

4 


Read  TITLE  nameliat 
Input 


Read  data  cards  to  gen- 
erate the  element  infor- 
mation 


Read  data  carde  to  generate 
the  collocation  points,  pivotal 
element  number  and  the 
mirror  image  t 


Cxhaost  N 
No ■ <^aU  influencing 
^^letneote  7 > 


Flf , 1 - How  Chart  for  Mata  Program 


(Continued) 


LENODE 

Array  containing  the  leading  edge  nodes 

LIST  1 

Input  namelist  name 

LMAX 

Number  of  entries  to  the  coefficient  table;  also 
the  number  of  nodal  information  data  card  to  be 
read  in 

MAXINT 

Number  of  integration  points  for  the  modified 
Gauss -Chebyshev  quadrature 

MIRKOR 

Node  number  of  the  mirror  image  point  with 
respect  to  the  centerline 

NEL 

Pivotal  element  number 

NELEM 

Number  of  elements 

NLE 

Number  of  leading  edge  nodes 

NMODE 

Number  of  mode  shapes 

NP 

Number  of  nodes 

QIMAG 

Array  of 

QREAL 

Array  of  QL 

SWPBK 

Sweepback  factor 

TITLE 

Array  containing  alphanumeric  information  for 
identification  purposes 

TOL 

Tolerance,  set  at  10  for  this  program 

UPWASH 

Array  of  the  upwash 

W 

Array  of  the  weights  of  the  modified  Gauss - 
Chebyshev  quadrature 

X 

Array  of  the  x-ordinates  of  the  nodes 

XEL 

Array  of  the  x-ordinates  of  the  elemental  nodes 

XK 

Reduced  frequency 

XLAMDA 

Sweptback  angle  of  the  leading  edge,  in  degrees 

XLIFT 

Array  of  the  lifts  at  the  nodes 

XM 

Mach  number 

XO 

Relative  position  in  x direction 

Y 

Array  of  the  y-ordinates  of  the  nodes 

YEL 

Array  of  the  y-ordinates  of  the  elemental  nodes 

YO 

Relative  position  in  y direction 

^ 


— 
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SECTION  IV 
INPUT  DESCRIPTION 


Input  cards  to  this  program  should  be  prepared  and  arranged  in  the 
order  described  below. 


A.  TITLE  CARD  (8A10) 


Information  on  planform,  Mach  number,  reduced 
frequency,  modes  and  mesh  spacing,  for  identifi- 
cation purposes 


B.  NAMELIST  Input 


$LIST  1 

XM 

Mach  number 

XK 

Reduced  frequency 

DELTA 

Mesh  spacing  as  measured  by  the  length  of  the 
side  of  the  characteristic  element 

MAXINT 

Number  of  integration  points  for  the  modified 
Gauss -Chebyshev  quadrature 

DEL 

Ratio  of  the  singular  strip  half  width  to  the 
element  half  width 

IMAX 

Maximum  number  of  characteristic  elements 
in  the  chordwise  direction  for  the  stencil 

JMAX 

Maximum  number  of  characteristic  elements 
in  the  Mach  line  direction  for  the  stencil 

XLAMDA 

Sweptback  angle  of  leading  edge  in  degrees 

NMODE 

Number  of  modes  in  the  set  of  mode  shapes 
(either  6 or  7) 

NP 

Number  of  nodes 

nelem 

Number  of  elements 

ICHECK 

Option  parameter  used  to  check  the  mesh 
correctness.  For  ICHECK  = 1,  a quick  run  is 
performed  to  print  out  the  element  information 
list,  as  well  as  a list  of  the  pollocation  points 
with  its  associated  influencing  elements.  For 
normal  run,  this  card  is  to  be  omitted. 
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NLE 

LENODE 

X 

Y 

$END 


Number  of  leading  edge  nodes 
Array  of  the  leading  edge  node  numbers 
Array  of  x-ordinates  of  the  nodes 
Array  of  y-ordinates  of  the  nodes 


C.  Card  for  Total  Number  of  Element  Information  Cards  to  Follow  (15) 

D.  Element  Information  Cards  (615) 

These  are  cards  to  generate  the  element  number  and  its  nodal  numbers 
in  a consecutive  manner.  Each  card  begins  a new  sequence. 


Col. 

I- 5 
6-10 

II- 15 
16-20 
21-25 
26-30 


Description 

Number  of  elements  to  be  generated  in  this  sequence 

Element  number  of  the  first  element  in  this  sequence 

First  nodal  number  of  the  first  element  in  this  sequence 

Second  nodal  number  of  the  first  element  in  this  sequence 

Third  nodal  number  of  the  first  element  in  this  sequence 

Fourth  nodal  number  of  the  first  element  in  this  sequence 
(leave  blank  for  triangles). 


For  example,  the  card 


22  13 


12 


generates  the  following  information 


Element  No. 


Node  Numbers 


9 

10 

11 

12 

13 


22 

13 

6 

12 

23 

14 

7 

13 

24 

15 

8 

14 

25 

16 

9 

15 

26 

17 

10 

16 

9 


while  the  card 


1 14  27  18  17 

generates  the  information  on  a single  triangle 


Element  No. 
14 


Node  Numbers 
27  18  17 


Element  nodes  are  ordered  in  a counterclockwise  direction,  starting  with  the 
most  downstream  node. 


E.  Card  for  Total  Number  of  Node  Information  Cards  to  Follow  (15) 


F.  Collocation  Point  Information  Cards  (415) 

These  are  cards  to  generate  the  collocation  point,  pivotal  element  number 
and  mirror  image  node  number  in  a consecutive  manner.  Each  card  begin,  a 
new  sequence. 


Col, 

1-5 

6-10 

10-15 

16-20 


Description 

Node  number  of  the  first  collocation  point  in  this  sequence 
Pivotal  element  number  containing  the  collocation  point 
Node  number  of  the  mirror  image  point 
Number  of  collocation  points  to  be  generated 


For  example,  the  card 


49  30  49  6 

generates  the  following  information. 


Node 

Pivot  Element 

Mirror 

49 

30 

49 

50 

31 

48 

51 

32 

47 

52 

33 

46 

53 

54 

34 

35 

45 

44 

10 
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while  the  card 


36  20  35  5 


generates  the  following 


Node 

Pivot  Element 

Mirror 

36 

20 

35 

37 

21 

34 

38 

22 

33 

39 

23 

32 

40 

24 

31 

SECTION  V 
SAMPLE  RUN 

The  case  of  a rectangular  planform  with  A = 2,  M = 1.2  and  k=  0.3  is  used 
to  illustrate  the  use  of  this  program.  The  planform  with  the  node  numbers 
and  element  numbers  is  set  up  as  in  Fig.  2.  For  a production  run,  a much  finer 
mesh  can  be  used,  and  the  main  program  can  be  dimensioned  accordingly.  The 
input  deck  for  this  problem  is  listed  on  page  17.  Some  suggested  values  for 
the  parameters  are:  MAXINT  > 12  and  DEL  = 0.7.  These  were  determined 
through  an  accuracy  study. 
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- 1 . . - ,66666666b 7 . - . 333 3333333 . 0 . . . 3 333333 333 » *b*b6g3 3333*3*3 , j . 

- 1 . . - . 83333  333  3 3.  - . 5 . - . 1 666666667  . . 1 666666667  . . 5 . . 8333333333 


SEND 
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2 
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13 

20 
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13 
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21 

15 

8 

14 
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20 
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14 

5 

21 
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22 

15 

21 

1 

26 

33 
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26 

6 

27 

34 

28 

21 

£7 

1 

33 

40 

34 

27 

5 

34 

41 

35 

28 

34 

1 

39 

46 

33 

39 

6 

40 

47 

41 

34 

40 

1 

46 

53 

47 

40 

5 

47 

54 

48 

41 

47 

1 

52 

59 

’♦6 

52 

6 

53 

61 

54 

47 

53 

7 

59 

60 

61 

53 

9 

1 1 

4 

10 

3 

17 

10 

17 

3 

24 

17 

23 

3 

30 

23 

30 

3 

37 

30 

36 

3 

43 

36 

43 

3 

50 

43 

49 

3 

56 

49 

56 

3 

64 

56 

63 

3 
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For  the  initial  run,  it  is  desirable  to  check  the  correctness  of  the 
mesh  first.  For  this  purpose,  a card  containing  ICHECK  = 1 can  be  included 
within  the  input  namelist.  This  option  will  bypass  all  the  time  consuming 
computations  and  output  only  the  mesh  information.  The  user  can  compare 
this  information  against  the  mesh  pattern  to  eliminate  any  input  error.  The 
output  for  ICHECK  = 1 for  this  sample  run  is  listed  in  the  following  pages. 


RECTANGULAR  A*2  • i M«l*2.  K*«3.  MODES*  1»  *•  *#x»  V*Y* 


X«K  »Y»T  » V i XT 


MESH  CHECK 


NODE 


X-ORD I NATE  T-ORDINATE 


1 

2 

3 
9 
S 

4 

7 

8 
9 

10 

U 

12 

13 
19 
15 

14 
17 
!• 

19 

20 
21 
22 

23 

29 
25 

24 

27 

28 

2» 

30 

31 

32 

33 
39 
35 

34 

37 

38 

39 

90 

91 

92 


• *50000*00 
•*50000*00 
••50000*00 
• *50000*00 
••50000*00 
• *50000*00 
-•50000*00 
••38996*00 
••38995*00 
••38996*00 
••38996*00 
••38995*00 
••38995*00 
••27889*00 
-•27889*00 
-•27889*00 
-.27889*00 
••27889*00 
••27889*00 
.•27869*00 
•• 14839*00 
••14839*00 
••14839*00 
••14839*00 
-•14839*00 
••14839*00 
••57763*01 
•.67783*01 
-.57783*01 
-•57763-01 
57783-01 
••57783*01 
••57783*01 
•52771*01 
•52771*01 
.52771*01 
•52771-01 
.Si*77  1-01 
•52771-01 

• 14332*00 

• 14332*00 

• 14332*00 


-•10000*01 

••44447*00 

••33333*00 

•00000 

.33333*00 

•44447*00 

•10000*01 

••83333*00 

-•50000*00 

••14447*00 

•14447*00 

•50000*00 

.83333*00 

-•10000*01 

••44447*00 

••33333*00 

•00000 

•33333*00 

•44447*00 

•10000*01 

-•83333*00 

-•50000*00 

-•14447*00 

• 14447*00 
•50000*00 
•63333*00 

-•10000*01 

••44647*00 

-•33333*00 

•ooooo 

•33333*00 

•64447*09 

• 10000*01 
-.83333*00 
-•50000*00 
••16667*00 

•16667*00 

•50000*00 

•83333*00 

••10000*01 

-•66447*00 

-•33333*00 
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93 

•14332*00 

•00000 

HH 

• 14332*00 

•33333*00 

46 

•16332*00 

•44447*00 

44 

•14332*00 

• 10000*01 

•♦7 

•27388*00 

••83333*00 

48 

•27388*00 

••60000*00 

49 

•27388*00 

••14447*00 

60 

•27388*00 

• 14447*00 

S| 

•27388*00 

•50000*00 

52 

•27388*00 

•83333*00 

S3 

• 38993*00 

•• 10000*01 

54 

•38993*00 

••64467*00 

55 

•38993*00 

••33333*00 

54 

•38993*00 

•00000 

57 

•38993*00 

•33333*00 

68 

• 38*493*00 

•64467*00 

69 

•38993*00 

• 10000*01 

60 

•99999*00 

••10000*01 

4 1 

•99999*00 

••83333*00 

62 

•99999*00 

••60000*00 

43 

,99999*00 

••14667*00 

44 

•99999*00 

•14467*00 

46 

•99999*00 

•60000*00 

44 

•99999*00 

•63333*00 

47 

•49999*00 

•10000*01 

6LEHCNT 

1 

2 

NODES 

8 2 1 

0 

9 3 2 

0 

3 

10  9 3 

0 

4 

11  5 *♦ 

0 

6 

12  6 6 

0 

4 

13  7 6 

0 

7 

19  8 l 

0 

8 

15  9 2 

6 

9 

14  10  3 

9 

10 

17  11  •» 

10 

1 1 

18  12  5 

1 1 

12 

19  13  6 

12 

13 

20  7 13 

0 

19 

21  18  * 

14 

16 

22  16  9 

15 

14 

23  17  10 

1 6 

17 

29  18  11 

17 

18 

25  19  12 

18 

19 

26  20  13 

19 

20 

27  21  1H 

0 

21 

28  22  15 

2 I 

22 

29  23  16 

22 

23 

30  29  17 

23 

2M 

31  25  IB 

29 

17 


25 

24 

27 

28 

29 

30 

31 

32 

33 
3*1 
35 

34 

37 

38 

39 
HO 
HI 
H2 
Hi 
HH 
H5 
H4 
H7 
H8 
H9 

50 

51 

52 

53 
5H 
55 

54 

57 

58 

59 

40 

41 

42 

43 
4H 
45 


NODE 


11 

12 

13 


32 

26 

1 V 

O 

33 

20 

26 

0 

3H 

28 

21 

27 

35 

29 

22 

28 

34 

30 

23 

29 

37 

31 

2H 

30 

38 

32 

25 

31 

39 

33 

26 

32 

HO 

3H 

27 

0 

HI 

35 

28 

3H 

H2 

36 

29 

35 

H3 

37 

30 

36 

HH 

38 

31 

37 

H5 

39 

32 

38 

H4 

33 

39 

0 

H 7 

HI 

3H 

HO 

H8 

H2 

35 

HI 

H9 

H3 

36 

H2 

50 

HH 

37 

H3 

51 

H5 

38 

HH 

52 

H6 

39 

HS 

53 

H7 

HO 

0 

5H 

H8 

HI 

H7 

55 

H9 

H2 

HB 

54 

50 

H3 

H9 

57 

51 

HH 

50 

58 

52 

H5 

51 

59 

H 6 

52 

0 

41 

SH 

H7 

53 

42 

55 

H8 

5H 

43 

56 

H9 

55 

4H 

57 

50 

56 

45 

58 

51 

57 

44 

59 

52 

58 

40 

61 

S3 

0 

41 

42 

5H 

0 

42 

43 

55 

0 

43 

4H 

54 

0 

4H 

45 

57 

0 

45 

44 

58 

0 

44 

67 

59 

0 

NEL  MIRROR  INFLUENCING  El CMENTS 


5 


Si*. 


After  the  mesh  has  been  verified,  the  program  can  be  run  without  the 

ICHECK  ~ 1 option  card.  The  output  consists  of  the  arrays  of  complex  lifts 

due  to  a unit  displacement  in  the  respective  mode  shapes.  The  output  is 

formatted  such  that  each  line  contains  four  nodal  lift  values,  with  the  real 

and  imaginary  part  given  in  pairs  (pages  24  through  27).  Following  these 

are  the  tables  of  generalized  force  coefficients  (page  28).  The  upper  table 

represents  Qi.  while  the  lower  table  represents  Q'.'..  The  output  for  the 
ij 

sample  run  is  listed  in  the  following  pages. 
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Appendix  A 


This  appendix  contains  a brief  outline  on  the  purpose,  method  and  use 
of  each  of  the  eight  subroutines.  The  principal  input  and  output  variables 
are  described.  The  subroutines  are  arranged  in  alphabetical  order  as  follows 

KERNEL 
LG SPAN 
POLYGN 
QIJ 

SGRHBS 

SGTRGL 

SINGUL 

TABLE 


SUBROUTINE  KERNEL 


PURPOSE: 

METHOD: 


USE: 


y 


2 

To  evaluate  the  reduced  kernel  function  K = yQ  K/2, 
where  K is  the  supersonic  kernel  function. 

The  nonplanar  form  of  the  oscillatory  sup  or  sonic 
kernel  function  was  derived  by  Harder  and  Rodden 
and  was  reduced  to  the  planar  form  by  A.M. 
Cunningham  in  the  appendix,  J.  Aircraft,  Vol.  11, 

No.  10,  October  1974,  pp.615. 

CALL  KERNEL  (XO.-YO,  XKREAL,  XKIMAG) 

Input 

XO,  YO  Coordinates  of  the  point  of  evaluation 
relative  to  the  collocation  point. 


Output; 

XKREAL,  Real  and  imaginary  parts  of  the  reduced 
XKIMAG  kernel  function. 


f*  . IP1PU.  J-:,  UWIPI'P  „ niyu 


SUBROUTINE  LGSPAN 


PURPOSE: 


METHOD: 


To  integrate  a function  with  inverse  square  singularity 
over  half  of  the  singular  strip  by  extracting  the  Cauchy's 
principal  value. 


Since  a singular  element  contains  either  the  complete 
strip  or  half  of  the  strip,  it  is  more  convenient  to  treat 
only  half  a strip  at  a time.  As  the  inverse  square 
singularity  occurs  only  in  the  spanwise  integration, 
the  chordwise  integration  can  be  performed  first  as 


ib(r>) 

F(n>  = / f(£,  r,)  de 

Sa(n) 


£b(r>)  • 5a(n)  V 


£wi  fi 


with 

and 


= f(Si(  n) 

A ^b(ri)  " Sa(ri)  . ^b(r})  + 

*i  = 2 Ci+ 2 


where  w.  and  are  some  Gaussian  weights  and 


abscissas  over  (-1,  1), 


A sixth  degree  quadrature  based  on  Lagrangian  inter- 
polation in  conjunction  with  the  Cauchy's  principal 
value  was  devised  by  Watkins  as 


r 

v-  r 


<y-<\) 


100c 


[l3(F1  + F?)  + 72(F  ^ +F6) 


+ 495(F3  +F5)  + (-1360)  F4] 


For  use  in  this  subroutine,  the  above  Quadrature  is 
modified  to 


29 


fr.,,  n 


USE: 


SUBROUTINES 

CALLED: 

ERROR 

RETURNS: 


/ 


llnl 


2 'ly  = 100 


k [nFi 


(y  - n) 

y-e 

for  the  left  half,  and  to 

y+e 


+72F2+496F3 - 680F4 


/ 


-Lin), 
(y-  n)‘ 


dy  = 


— - — Tl3F 
100  e L 7 


+ 72F6  +495F5  - 680F4] 


for  the  right  half. 

CALL  LGSPAN  (X,  Y,  XEL,  YEL,  EINT,  W.MAXINT,  F) 
Input; 

Coordinates  of  the  collocation  point. 


X,Y 

XEL, 

YEL 

EINT 

W 


Arrays  of  nodal  coordinates  of  the  half 
strip. 

Array  of  Gauss -Chebyshev  quadrature 
abscissas . 

Array  of  Gauss-Chebyshev  quadrature 
weights . 


MAXINT  Number  of  integration  points. 

Output: 

F Array  of  integrated  values  at  nodal  points 

KERNEL 


None 
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SUBROUTINE  POLYGN 


PURPOSE: 

METHOD: 


USE: 


SUBROUTINE 

CALLED: 


To  perform  the  integration  over  a regular  triangular 
or  quadrilateral  element. 

The  arbitrary  triangular  or  quadrilateral  element  is 
mapped  into  a square  region  -1  < £ <1  and  - 1 < rj  < 1 . 
Integration  is  accomplished  by  a repeated  application 
of  the  Gauss -Chebyshev  quadrature  >n  both  directions. 
Since  the  mapping  is  different  for  a triangle  and  a 
quadrilateral,  different  computational  routines  are 
employed. 

CALL  POLYGN (X , Y,  XEL,  YEL,  EINT,  W,  MAXINT,  F, 
IANGLE) 


Input: 

X,  Y 

XEL,  YEL 
EINT 

W 

MAXINT 

IANGLE 

Output: 

F 


Coordinates  of  the  collocation  point. 

Arrays  of  nodal  coordinates  of  the  element. 

Array  of  Gaus s -Chebyshev  quadrature 
abscissas . 

Array  of  Gaus s -Chebyshev  quadrature 
weights. 

Number  of  integration  points. 

Number  of  angles  in  the  polygon  element. 


Array  of  integrated  values  at  nodal  points. 


KERNEL 


ERROR 

RETURNS:  None 
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SUBROUTINE  QIJ 


PURPOSE: 

METHOD: 


To  compute  the  generalized  force  coefficients  for  a 
given  load  distribution. 

With  the  finite  element  approximation,  the  general- 
ized force  coefficients  are 


Q..  = JJ L(x,  y)  n‘(x,  y)  dxdy  . \je) 


n=1  A<e> 


N 


where  denotes  summation  over  all  elements, 


n = l 


f are  the  modes,  N*  are  the  shape  functions  and 
i 

X<e)  is  the  elemental  lift  vecvor.  The  real  and 
imaginary  parts  are  defined  as 


Qij  - Re<Qij> 


and 


Q|j  = Im(Q..)/k 

where  is  the  reduced  frequency.  The  integrations 
are  performed  using  the  Gaussian  quadrature  (2  points 
as  set  up  in  the  subroutine). 

CALL  QIJ(X,  Y,  INFO,  XLIFT,  FF,  QREAL,  QIMAG , 

NP,  NELEM,  NMODE,  MX) 


Input: 

NP  Number  of  nodal  points. 

NELEM  Number  of  elements. 

NMODE  Number  of  modes. 

MX  Maximum  number  of  nodec  the  the  elements. 

X,Y  Arrays  of  nodal  coordinates. 

INFO  Array  of  element  information. 

XLIFT  Array  of  modal  lift  distribution. 

FF-  Array  for  temporary  storage. 


Output: 


QREAL  Real  part  of  the  generalized  force 
coefficient,  QJ.. 

QIMAG  Imaginary  part  of  the  generalized  force 
coefficient,  QI'j* 


SUBROUTINES 

CALLED: 

None 

ERROR 

RETURNS: 

None 

PURPOSE: 


METHOD: 


USE: 


SUBROUTINES 

CALLED: 

error 

RETURNS: 


SUBROUTINE  SGRHBS 

To  perform  the  integration  of  a rhombic  element 
with  the  singular  strip  passing  through  either  the 
left,  middle  or  right  node. 

The  rhombic  element  is  divided  into  the  right  half 
triangle  and  the  left  half  triangle  . These  tria^la 
elements  can  be  either  regular  or  singular.  The 
integrations  are  accomplished  by  calling  another 
subroutine  and  the  results  are  re -assembled. 

CALL  SGRHBS(XNODE,  YNODE,  XEL,  YEL,  EINT, 
w MAXINT.F.DEL) 


Input; 

XNODE, 

YNODE 

XEL,  YEL 

EINT 

W 

MAXINT 

DEL 


Coordinates  of  the  collocation  point. 

Arrays  of  nodal  coordinates  of  the  element. 

Array  of  Gauss -Chebyshev  quadrature 
abscissas. 

Array  of  Gauss -Chebyshev  quadrature 
weights. 

Number  of  integration  points. 

Ratio  of  the  singular  strip  half  width  to 
the  element  half  width. 


Output: 

F Array  of  integrated  values  at  nodal  points. 


SINGUL 


None 


SUBROUTINE  SGTRGL 


PURPOSE: 


METHOD: 


USE: 


SUBROUTINES 

CALLED: 

ERROR 

RETURNS: 


To  integrate  a general  triangular  element  with  the 
singular  strip  passing  through  either  one  of  the 
three  vertices. 

A vertical  line  through  the  lower  vertex  divides  the 
element  into  two  sub -triangles,  which  can  be  either 
regular  or  singular.  The  integrations  are  accom- 
plished by  calling  another  subroutine  and  the  results 
are  re-assembled. 

CALL  SGTRGL(XNODE,  YNODE,  XEL,  YEL,  EINT, 

W,  MAXINT,  F,  DEL) 

Input : 

XNODE,  Coordinates  ol'  the  collocation  point. 

ynode’ 

XEL,  YEL  Arrays  of  nodal  coordinates  of  the  element. 

EINT  Array  of  Gauss -Chebyshev  quadrature 

abscissas. 

W Array  of  Gauss -Chebyshev  quadrature 

weights. 

MAXINT  Number  of  integration  points. 

DEL  Ratio  of  the  singular  strip  half  width  to 

the  element  half  width. 

Output: 

F Array  of  integrated  values  at  nodal  points. 

SINGUL 

None 
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PURPOSE: 


METHOD: 


USE: 


SUBROUTINE  SINGUL 

To  integrate  the  special  type  of  triangular  element 
bounded  by  a vertical  line  and  two  other  straight 
lines.  The  triangular  element  may  be  singular  or 
regular. 

The  triangular  element  is  first  tested  to  see  if  it  is 
singular  or  regular,  k or  singular  element  further 
test  is  conducted  to  locate  the  singularity  with  re- 
spect to  the  element  such  that  the  element  can  be 
divided  into  a singular  strip  and  a regular  polygon. 
Integrations  are  accomplished  by  calling  other  su 
routines  and  the  results  are  re-assembled. 

CALL  SINGUL{XNODE,  YNODE,  XEL,  YEL,  EINT, 
W,  MAXINT,  F , DEL) 


Input: 


XNODE, 

YNODE 

XEL,  YEL 

EINT 


MAXINT 

DEL 


Output: 


Coordinates  of  the  collocation  point. 

Arrays  of  nodal  coordinates  of  the  element. 

Array  of  Gauss-Chebyshev  quadrature 
abscissas. 

Array  of  Gauss-Chebyshev  quadrature 
weights . 

Number  of  integration  points. 

Ratio  of  the  singular  strip  half  width  to 
the  element  half  width. 


Array  of  integrated  values  at  nodal  points. 


SUBROUTINES 

CALLED: 


POLYGN, LG SPAN 


ERROR 

RETURNS: 


None 
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SUBROUTINE  TABLE 


PURPOSE; 


METHOD: 


USE: 


SUBROUTINES 

CALLED: 

ERROR 

RETURNS: 


To  create  a table  of  the  weighted  kernel  function 
coefficients  for  a given  uniform  characteristic 
mesh.  This  table  is  stored  for  later  table  look-up 
in  the  solution  process. 

A stencil  of  uniform  characteristic  mesh  large 
enough  to  cover  the  most  extreme  case  for  the 
plan  form  is  set  up.  Because  of  the  symmetry  in 
the  spanwise  direction,  only  the  rhombic  elements 
on  one  side  need  to  be  evaluated  and  stored.  Each 
element  is  uniquely  defined  by  a pair  of  relative 
indices  based  on  its  relative  location  from  the 
collocation  point. 

CALL  TABLE(COEF,  DELTA,  EINT,  W,  MAXINT, 
DEL,  IMAX,  JMAX,  LMAX) 


Input: 

DELTA 

Length  of  the  side  of  the  characteristic 
element. 

EINT 

Array  of  Gauss -Chebyshev  quadrature 
abscissas. 

W 

Array  of  Gauss-Chebyshev  quadrature 
weights . 

MAXINT 

Number  of  integration  points. 

DEL 

Ratio  of  the  singular  strip  half  width  to 
the  element  half  width. 

IMAX 

Number  of  elements  in  the  chordwise 
direction  for  the  stencil. 

JMAX 

Number  of  elements  in  the  Mach  line 
direction  for  the  stencil. 

LMAX 

Number  of  entries  to  the  table. 

Output: 

COEF 

Array  of  the  weighted  kernel  function 
coefficients 

SGRHBS,  POLYGN 


None 
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Appendix  B 

PROGRAM  LISTING  OF  FESKAP 


r>  r,  r . c a ('  n n r r.  n r r n 


wmm 


»uw— ■wmssm 


Copy  evailable  to  bdC  dees  not 
permit  fully  legible  reproduction 


Hr0GNAM  FeSKAFI  UPuT.GOTPuT,TAPeS«INPuT.TaPeo*uUTPUT) 

C R n.E*  FIMTi  li_E>.cnT  .nerNel  FUNCTION  .'H.TMWU  In  -URERsuNtC 
C LIFTING  bUHhAct  TntoKY 

C u-VeLOPeU  sY  V • Y»  C.  YDONG  «T  LOCKHEED  - HUNTSVILLE 
CgMPLEX  F » SUi’t  .UPWASH.CCEF  «XL  t F T 

DIMENSION  X<  67  I « Y < 67)  . INF0<65«4>  . Xl  IF  T ( 67*7 ) «COtF  < 4 .26  > • 

I LEN0Dt<7) . IbUF<65)  .FF<67.7> 

OIMeNoIU.m  F ( 4)  • Aee  < 4 ) « YtL  < 4 ) * GUM  17)*  UPS ASH<  7 ) » T 1 1 ll  ( b ) . 

1 oRtALf  7*7  » ,ulMAG<  7*7  > * w l 30 ) « 1 1 NT ( 30  > 

CgMnON  *■••»  XK  «3t  I As- .SET A 
egU  1 VALeNee  ( r r l i * 1 > * COeF  < I « I > > 

Trie.  FIRST  0 1 .'iti'YS  i 0,i|  oTATEMfcNT  Is  D l MENS  1 uNi-O  AG  X(NP>.Y(NF>. 
|.vrutNEUcl>.«4)  .XL  IF  M .r-.M-iOtje)  «COtF  U.u'iAXI  .LliY-UlIIVl  > . loU^tNELfc  > . 
Fr  (NP.NMODe) 

LRN  UL  ALTERED  BY  THE  USER  TO  FIT  THE  PkOcJLtF. 


IoUF  AND  FF  are  TtMPJRAHY  STOhAG^  ARRAYS 


GEL 

UluTA 
1 i-riH-N 


C l i-i  AX 
C 

C gr  AX 
L 

L LLNJUe 
C i'AXINT 
0 

c l.-iAx 
L WeLc'I 
C l«UL 
L NmOUc 


= WIDTH  OK  Trie  SINGULAR  STRIP  RELATIVE  TC  The  SIgTf,  CF 
The  ELFi-i-iMf 

* Length  jf  Trie  S I gl  OF  THE  Chmral Teh  1 ST  1 0 eeeMenT 
= DUMMY  R0n  Tg  CHECK  THe  MESH  eGRRECTNlSS 

IF  1 CHECK  IS  SeT  TO  I 
uTHeR^ISL  i CHECK  IS  SET  TC  0 

* MAXIMUM  rU  <iJER  OF  eHARACTer l sT 1L  ELtMeNT  IN  THe 
CHORD W I S E DIRECTION 

s MAXIMUM  nU.-ldeR  OF  CHARACTER  I sT  I o ELEmeNT  IN  THL 
MACH  LINE  DIRECTION 

r ARRAY  CDMaININO  THE  LE  AO  l NO  eUGE  NOOeS 
= NUMGeR  JF  1 ITeGkATION  POINTS  FOR  THE  MODIFIER 
GAUSS-CHeBToHeV  DoAJRATUHE  (MAXIMUM  SET  aT  JO  > 

= IMAX*JF4A-(  I MAX* < I MAX- I > >/2 
= NO.  OF  ELEMENTS 
= NO.  Or  LlAdInG  eji)L  NOUES 
s NO.  OF  lOUt  SHAPeS 


C 


O 

C 


M-  = NO.  OF  RUOco 

XK  = REDUCED  Fi<lJUENLY 

xlamoa  s s^lptuaCk  angle  of  leading  edge  in  degrees 

XM  = MALH  NuMoLR 


data  tul/i.e-^/  . , 

NAMEL  1ST  /LlStl/XM.XK, DELTA. MAX I NT .DEL*  IrlAX.oMAX. XeAMJA .NwluDL *NP » 
1 rELeM « 1 CHECK . NeE . eeNODe . X . Y 

1 NDREL  (A.a)*lF"IXt  CST  l * ( A+bET  A*D  ) +.S  > ♦ l 
LINEAR!  I . 0 ) : ( ( l - 1 ) * < UMAX* UMAX— I > >/2+U 


c read  title  oaru  Aru  ramelist  input 

c 

HEAU(is»9200>  TITLE 
HEAD (5.LISTI ) 

DeTASg»XM#XM-1 . 

UETAsSOrKOeTASj) 

CST 1 =.S*XM/(BtTA»DELT  A 1 

IF  ( ICHeCa.EJ. 1 ) 4R  I T E < 6 * 96 1 0 ) TITLE 

IF  ( ICHeCK.EO.1  ) a'RITE(€.W210)  C I . X ( I)  . Y ( l ) . I = l «NP  > 


A Do  I 

A WLC 

A OuJ 
A OyA 
A OCS 
A 00b 
A 007 
A 00B 
A 009 
A 010 
A 0 1 1 
A OIL 
A y i j 
A 01 A 
A vis 
A 0 1F 
A 017 
A Slo 
A 0 19 
A CcO 
A CL  I 
A CeC 
A OeJ 
A OeA 
A 0eo 
A Ocu 
A 027 
A Ocu 
A OtV 
A 0g0 
A 0-;l 
A 0j2 
A OeD 
A OgA 
A OJa 
A 0gO 
A 0g7 
A 0 go 
A C g9 
A 0 AO 
A 04  1 
A CAL 
A 043 
A 044 
A 0 AS 
A 046 
A 047 
A 040 
A 04  9 
A DSC 

F4  w w A 

A 

A wE3 
A OCA 
A 0‘g3 
A 056 
A C‘j7 
A OS  (J 


C GlNCRATE  fc.Lt.MENT  INFORMATION  FROM  OATA  CARDS 
C 

READ ( 5 • 9000 ) UMAX 

IF  ( ICHECK*EO«l t WRITE (6*9620) 

DO  12  K«1*KmAX 

READ) 5*9000)  IRtPT* I 1 *N1 *N2*N3*N4 
00  10  l-LIREPT 
IMl-l-l 
12-1 1 ♦ I Ml 

INFO!  I2*1)>n1+1M1 
INFO) I2*2)-N2+IM1 
INFO!  I2*3)-n3+1M1 
INFO! I2*4)«N4-MM1 
IF  (N4-CQ-0)  1NF0< I2«4)-0 

IF  I ICHuCK.EG-l)  WR I TEI 6*9600 ) 12 • INFO( 12 • 1 ) * INFO! 1 2 *2 ) . 

I INFO ( 12*31*1 NFO( 12  *4 ) 

10  CONTINUE 
12  CONTINUE 

IF  ( I CHECK .EO* 1 ) oo  TO  5001 

c 

C OtNERATE  THc.  MOD IF lc.3  GAUSS-CHEBYSHEV  WEIGHTS  AND  ABSCISSAS 
C 

CST-3.14159265359/MAXINT 
DO  20  I-1*MAX1NT 
ARg*< l-*5)*CST 
EInT( I l-COS(ARG) 

20  Wl I >»CST*SIn< ARG) 

c 

C CLEAR  the  NODAL  LIFT  VALUES 
C 

DO  30  MODE-l »NMODE 
DO  30  I - l • Np 

30  XLIFT(I*M00e)-CMPLX(0«*0.) 

c 

c generate  the  weighted  kernel  function  coefficients  for  the  stencil 

C OF  UNIFORM  CHARACTERISTIC  MESH 

c 

UMAX -LI NEAR ( I MAX (UMAX) 

CALL  TABLE ( COfcF • DELTA  *E I NT • W * MAX  I NT *DEL  * IMAX • JMAX • LMAX ) 

C 

C COMPUTE  THE  SWEEP6ACK  FACTOR 

C FOR  SUBSONIC  LEADING  EDGE*  SET  SWEEPBACK  FACTOR  TO  1* 

C 

CsT- TAN( • 1 745 3292BE- 1 *XLAMOA ) /BETA 
SWPBK-l . 

CST* 1 »-cst*cst 

IF  (CST.GT. 1 »t-10)  SWPBK- 1 »/SURT  ( CST ) 

CST--2**SWPbX/8ETA 

c 

C COMPUTE  THE  STARTING  SOLUTION  AT  THE  LEADING  EDGE 
C 

00  40  I-1*NlE 

1 1 -LENOOE ( I ) 

XNOOE-XI I 1 ) 

YNOOE-V« I I ) 

XlSQ-XNOOE*XNODt 

VlSO»YNODE*VNUDc 


059 

060 

061 

062 

063 

064 

065 

066 

067 

068 

069 

070 

071 

072 

073 

074 

075 

076 

A.. 077 

078 

079 

OtoO 

081 

082 

083 

084 

085 

086 

087 

088 

089 

090 

091 

092 

093 

094 

095 

096 

097 

098 

099 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1 10 

1 1 1 

1 12 

1 13 

1 14 

1 15 

1 16 

j 


i 


< 


i 


40 


| 

4 


d 


non  n n n non 


XL  I F T < I 1 « 1 ) >C3T * CMPLX < 0 • * XK l 

XL  IF  T I 1 I *2)  *CST*CMPLX(  1 « *X<*XNODE  ) 

XL IFTI I 1 . 3>*CST*CMPLXIXN0UE*XN0DE.XK*X1SU) 

XL l F T I I 1 *4) «CST*CMPLX(0**XK»YlSU) 

XL l F T < I 1 *5)  =CST*CMPLX<2.*XNODE*YISO*XK»XISU*YISO) 

XL  IFTI I 1 • NMODt— 1 ) “CST *CMPLX I 0»  * XK# VNODE  > 

XLIFT<  I 1 iNMOJl  ) *CST*CMPLX<  YNOOt  * XK*XNOOc* YNOOE ) 

40  CONTINUE 

5001  CONTINUE 

c 

c loop  Through  the  unknown  nooal.  points  to  compute  the  lift 
c 

«EAD (5*0000)  LMAX 

IF  < ICHECK«E(i«l ) WR 1 TE (6*99 10) 

00  1000  LN00E“1*LMAX 

HEAD < 5*9000 ) N1«N2*N3. INEPT 
DO  1000  1*1. I kept 

1 MI  * I — 1 
NODE*Nl+IMl 
NEL  * N2  ♦ I M 1 
MIPH0H-N3-1M1 
XNOOE  * X ( NODE  I 
YNOOE* V (NODE ) 

IF  ( ICHECK.ECI«1 ) GO  TO  5002 
C COMPUTE  THE  MOOAL  UP WASH  AT  THE  COLLOCATION  POINT 
XI SU*XNOOE*XNODt 
Y I SO* YNOOE* Y NODE 
UPWASHI I >*CmPLXIO»*XK) 

UPW ASH ( 2 ) * CMP LX (1 • • XK*XNODE ) 

UPW ASH 1 3 ) * CmPLX I XNOOE +XNODE • XK*  X I SO ) 
UPwASH(4)*CmPLXI0«*XK*YISU> 

UPW A SH ( 5 ) * CMP L X ( 2 . * XNODE*  Y I SO  * XK*X  I SO*  Y I SO ) 

UPWASHI NMOOE* 1 1 * CMPLX ( 0 • • XK*  YNOOE ) 

UPWASHI NMOOE  ) “CMPLX I YNOOE*  XK«XNOOE* YNOOE  > 

CLEAN  THE  SUMS  BEFORE  ACCUMULATION 

OO  100  MQOE*l • NMOOE 
SUM! MODE )*CmPLX(0* *0* ) 

100  CONTINUE 

5002  CONTINUE 

LOOP  THROUGH  THE  FORWARD  ELEMENTS 
KOUNT  *0 

00  BOO  L* 1 *NCL 

1 angle *4 

IF  I INFO(L*4)*EU*0)  1 ANGLE *3 

ASSIGN  THE  ELEMENT  NOOAL  POSITION 

OO  200  K* I « I ANGLE 
KK» iNFOILtK) 

XEL(K)*X(KK) 

200  YEL(K)*Y(KK) 

xo*xnode-xel ( 1 ) 
yo*ynooe-yel ( 1 1 
c 


A 117 
A 1 IB 
A 1 19 
A 120 
A 121 
A 122 
A 123 
A 124 
A 125 
A 126 
A 127 
A 126 
A 129 
A 130 
A 131 
A 132 
A 133 
A 134 
A 135 
A 136 
A 137 
A 1 38 
A 139 
A 140 
A 141 
A 142 
A 143 
A 144 
A 145 
A 146 
A 147 
A 148 
A 149 
A 150 
A 151 
A 152 
A 153 
A 154 
A 155 
A 156 
A 157 
A 158 
A 159 
A 160 
A 161 
A 162 
A 163 
A 164 
A 165 
A 166 
A 167 
A 168 
A 169 
A 170 
A 171 
A 172 
A 173 
A 174 


I 


5 

I 


9 


! 


i 

i 


I 

I 


i 
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C SKIP  THfc  1LST  Fed  THt  PIVOTAL  fcLEMtNT 
C 

IF  4L.LO.NtL)  GO  TO  201 

C SKIP  THt  ELtMENT  IF  IT  IS  NOT  WlTHlN  THE  FORWARD  MACH  CONE 

IF  4 4XO-BETA*AB34YOl I .LT.4-T0L)  I GO  TO  BOO 
201  CONT INUE 

K0UNT.K0UNT-.1 
lOUFIKOUNT ) *L 

IF  4 ICHtCK.EO.l I GO  TO  BOO 
IF  4 1 ANGLE . EO . 3 ) GO  TO  400 
C 

' WtGULAR  CHARACTER  1ST ic  ELLMENTS 
C COMPOTE  THE  HELATIVt  INDICES 
C 

LL.INDREL4XO.YO) 

MM.1NDREL4XO.-Y0) 

IF  4LL.GT.MM)  GO  TO  300 
C 

C LOWER  TRIANGLE  OF  TABLE 
C 

Ll.LINEAH4LL.MM) 

F4ll.COEF41.Ll) 

F42I.COEF42.L1) 

FOI.COEF43.L1) 

F44l.COEF44.Ll) 

GO  TO  600 
SCO  CONTINUE 

c 

C UPPER  TRIANGLE  OF  TABLE 
C 

L1.L1NEAH4MM.LL) 

F C 1 1 >COEF(  l.Ll) 

F42l.COEF44.Ll) 

F40l.COEF4J.Ll  ) 

F 4 4 ) "COEF 4 2. L 1 ) 

GO  TO  600 
C 

c TRIANGULAR  FILL-IN  ELEMENTS 
C 

400  CONTINUE 
C 

C TEST  FOR  SINGULARITY 
C 

DO  BOO  K* 1 . 1 ANGLE 

IF  4 ABS4  YNOo£-YtL4KI  ) .LT.TOL)  GO  TO  S02 
BOO  CONTINUE 
C 

C NON-SINGULAR  POLYGON 

CALL  P0LYGN4  XNOilE*  YNOOE  . XEL  * YEL  «E  INT  «W  .MAXINT  .F  • I ANGLE  ) 
GO  TO  600 
C 

C SINGULAR  TRIANGLE 

C 502  CALL  SGTRGL4XNODE«YNOOE. XEL. YEL.EINT.W. MAXINT. F .DEL) 

600  CONTINUE 


A 175 
A 176 
A 177 
A 17B 
A 17V 
A 1B0 
A 1B1 
A 1B2 
A 1B3 
A 1B4 
A 1B5 
A 1B6 
A 167 
A 166 
A 169 
A 190 
A 191 
A 192 
A 193 
A 194 
A 195 
A 196 
A 197 
A 19B 
A 199 
A 200 
A 201 
A 202 
A 203 
A 204 
A 205 
A 206 
A 207 
A 20B 
A 209 
A 210 
A 211 
A 212 
A 213 
A 214 


A 215 
A 216 
A 217 
A 216 
A 219 
A 220 
A 221 
A 222 
A 223 
A 224 
A 225 
A 22 6 
A 227 
A 226 
A 229 
A 230 
A 231 
A 232 
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o r.  n 


wmgmzmw 


Copy  available  to  DDC  does  not 
permit  fully  legible  reproduction 


COPY  AVAILABLE  TO  DDC  DOES  NOT 
...  ,00  ... .,o..ooC  PERMIT  FULLY  LEGIBLE  PRODUCTION 

^cINKOlL.M  " " 

DO  700  WOUtr  1 »N.-, JOe 

9UMl  ,-iODt  I = sUi'i  ( MODE  J + XLIF  T (<K  « F,UDE  > *F  IK  > 

70  C CONTINUE 
dOCl  Continue 

lr  ( 1 Chu wK • \L  • 1 ) oU  TO  eOQ3 
WM|  Tt<6*y<»20)  NUjc  « NtL  «l"i  1 KWOK 
MK| TE I 6.0000  > t I OUF (<  ) .K* I .KOuNT  > 

Go  To  luoO 
l,CUJ  ConT  INUE 

O ruCAo->E;  T He.  HIVoTmL  r.  eLMENT  *AS  T«fc.ATtD  UAoT 

C r ( 1 ) CONTAINS  T Ht  COEFFICIENT  FOR  TMl  UNKNOWN  P I VOTAo  LIFT 
0 joObT  i TUTL  in  V- 1 r-  T r>  wM.O,  InAOe  TakINo  AuVmnTmoe  OF  SY.W-TA  Y oW 
C A^TI-SY/'IETHY 

c 

DO  -J 00  NODE  r 1 » iNi'l  JOL 

XulFT(Nu0E.VoUi)*l0.2«31Ob3071b*UPWAi,H(,v,00t  J-SUi'IWOt)  > /r  I 1 ) 

Xu|rT  (Ml  N«OR  ..'iJUt  > 1 XL  1FT  (NODt.F.UoC  ) 

•JuC  Continue 

Xl  | F T ( M I .«« oR  • . .MODE - I ) =-Xe  1 F T ( >•« I K«J»  • Nl»,wUE - I > 

Xu  I F T (M|r*<Ui?t..riUDL  >*-XL1FT(n1nNOH.NmuUE ) 

1000  LONTINUL 

Ir  (ICHt-CK.Fo.il  oTuP 
DO  2000  ,*IOoF  = i • ,VIUUL 
Wh  ITL(6«Ol  00  I FiOJe 

WHITE  (6.05C0)  ( XLlFT(  | . f*,oufc. I « l* l «NP> 

2GoO  LOnTINOL 

ooi*iPwTe  T He.  i luw  F oHuu  lOeFf  IuIliyTj*  -.Id 


call  uI0( X.Y.  t HF  0 * XL  I F T « f F « OKu  AL  « o I fiAG  • Ni-  » NuuC  •• « m, 'loot.  A ) 

W«lTc(6.0300)  TITLE 

IF"  (NMOOe.EO.O)  WR1 TE (6.OO10 > 

IK  ( Nw.oDt.LO.  7 I whI  TE  (tj.voOO  I 
Do  3000  1*1.  Ni'iU  Jt 

WHI  TE  ( 6.9A0C  I I » ( UWEALI  I .0)  .0*  1 .NfioDE  ) 

3000  CONTINUE 

«r|  fc( 6*0700  I 

IF  (NMOOe»cQ«6)  wN  I Tt  ( 6*06  1 0 ) 

IF  ( NMOOc.CG.  7 I WR l Tt< 6*0600) 

DO  3 100  1*1.  Ni'iJOfc 

W«  I TE  ( 6.0 A Of!  I 1 . ( L iF.AGt  I .0  I . 0*  l .NMOOE  I 
3100  CONTINUE 
oToP 

0000  FORMAT  (61b) 

0100  FORMAT  ( IHI  . 34nCo,-.PLEX  L IKT  DISTHIuUTlON  FUh  f.vu..  . I 2 . / > 

*2uU  FOHNAT(dAIO) 

Y2I0  FORMAT! IX. l5.10X.2Elb.31 
YoOC  KORfiAT  ( I H 1 . ///  . jX«OA10//''//> 

»AoO  FORMAT!  1H0*  I 5 » 1f*7e  15.A  ) 

ybOO  FORMAT <2E  1 b«  6 * AX  ict  1 b .6*  AX  « 2E  Ib.d.AX.  iib.  I b.B  ) 

9bo0  F'URMATI  IX.  15. 10X.AI5) 

YbIC  FuRmATI  IHI  ./✓/ • iX.oAlO.////.  IX.  ISHMtbH  //✓//. ~X» 

I AHNODfc  » 1 ?X  . 1 Of(X-OkO  INATc.  . bX  . 1 UHY-UHu  I NaTe»  //  • 

VOEO  FORMAT!/////. IX.  THtLfcMENT.VX.bHNOUtis.// > 

0700  FORMAT!////) 


A c J J 
A 22A 
A 235 
A 

A 237 
A 23U 
A 23/ 

A 2 AO 
A CA  1 
A 2A  2 
A 2A3 
M c A A 
A 2Ab 
A c AC 
A u A 7 
A uAo 
A _ A j 
A clC 
A w. i 
A ? - ? 

A 253 
A 2 -A 
A f..o 
A 2 30 
A CL  7 
A l-C 
A 2-  J 
A 200 
A col 
A cu2 
A CO  J 
A 2oA 

A co6 
A co7 
A coo 
A cC  f 
A 270 
A c7  1 
A 27c 
A 273 
A 2 7A 
A c 75 
A 276 
« 277 
A c7o 
u c7v 
A c3w 
A co  1 
- coc 
A c 03 
A cdA 
A 2d5 
A 2o6 
A co7 
A cob 
A „o  > 
A 
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9«OU  FORMAT.bAMHl.BA.bHJ  ■ | . lOK.bHJ  - 

i-w.cMi  m 5i10K«3HJ  * 6«  lOXtbHJ  * rf  1 
98 1 0 FORMATOXMHl.8X.bHJ  - I . IQX.bHJ  - it.IOX.bHJ  - JMOX.bHJ 

l lOXtiHJ  * b.tOXOHJ  « 6/) 

9900  FORMAT. 32X .25 14)  iciluU.bX  • 

99 | o FORMAT. ////////. 3X.4HN00E.6X.3HNfcL.AX.6HMlfi«0H.bX. 

I 20H|NFLUlNCING  tutMCNTb*// > 

99<iO  FORMAT. 1X« 15*21 10) 


3UBK0UT 1 Nfc  K. £ KNtL  t XO  • YO  • XKWfc  Ai_ • XX 1 MAG ) 


C SUPERSONIC  KERNEL  FUNCTION  OENlVEU  bY  HARDt-H  ANJ  WOUOcN 
C AS  GIVEN  BY  A M CUNNINGHAM  IN  APPENDIX  OF  J*  AIRCRAFT. 

C VOL*  11*  NO*  10*  197A* 

C 

DIMENSION  Allll 


COMMON  XM«  XK  *dfc T ASG • BETA 

— — OAT* Wt^<24-lJl6l9b«t*796d027.-E4*991079*  111  • b*» I V6 • 

1 UOa*  73XBB  » A 1 • 1^630^-545.  *9d537*  644  *761  ad.  “ 

d 64. 27901  j./  ~~~~~ _ 

DATA  C/.372/ 

R-SURTIXO*XO"btTASU*YO*YO> 

IF  (XK.LT.l *E-Dl  00  TO  400 

AVO-ABS< YOI 

ARG-XK*XO 

CS-COSt ARG) 

bNa  S I N ( AhIG  ) 

IF  < AYO*LT • 1 • t-& ) GO  TO  SOO 


-271*43549. 
326* 7275b. 


C Gc.RN£HAL  FURM  OF  NtKNEL  FUNCTION 
C 

XKAYO«X**AYO 
XKY0SU>XKAY0*AAAY0 
b2Y0  1 ■ 1 ./«  BETASU*AYO  I 
XMR-XM*H 

U1«B2Y01*«X0“XMKI 
U2-62Y0 1 * < X0+ XMK I 
£ 1 -EXP  t -C*AbS ( U1 1 I 
E2-EXPt-C*U2 I 
CST-O* 

CST 1 ■ I • 

CST2-1* 

SUM-0* 

SUM1RL-0* 

SUM  I IM-O* 

SUM2RL-0. 

SUM2IM-0* 

DO  100  1 « 1 * 1 1 
CST-CST+C 
CST l -CST l *E | 

CST  2-CsT2*E2 

COEF- At  1 1/ 1 CST*CST  4XKY0SU ) 

SUM - SUM+COEF 
CO£FRL-CST*C0tF 
COEF  IM— XXAY0*CUEF 
SUMIHL«SUMIRL+CUEKWL*CST l 
SUM  1 IM-SUM1 JM+COfcF IM*CST 1 
SUM2HL-SUM2Ru+CUtFWL*CST  2 
SUM2 IM-SUM2 l M+COEF IM*CST2 

|U0  continue 

A*G-XKAYO*Ul 
CS| -COSt  ARG) 

SNt-SIN(ARG) 

AHG-XKAY0*U2 
CS2-C0St ARG) 

SN2-S INt  ARG) 

X 1 1 2HL-XKAYQ* I SN2#SUM2RL-CS2*SUM2  1 M ) 


a 001 
3 002 
3 003 
3 004 
3 005 
3 006 
i 007 

a oob 

a 009 
a CIO 

3 0)1 

- ou 

3 013 
o C 1 4 
3 CIS 
s 016 
3 017 
u 016 
d 019 
u 020 

d 021 

a 022 
a 023 
3 024 
u 025 
a 026 
d 027 
6 026 
a 029 
b 030 
w 031 
a 032 
o 033 
a 034 
a 035 
a 036 
3 037 
o 036 
6 039 
a 040 
3 04  1 
a 042 
6 043 
b 044 
b 045 
a 046 
d 047 
d 046 
d 049 
b 050 
b 051 
b 052 
6 053 
b 054 
d 055 
b 056 
3 057 
d 056 


45 


non  r r *» 


xl  1 21  M*XXAYO*  ( C«2*50M2RC*aN2*i»UM2  IM  > 

CST l*SNl*6UMlRL 
CST2*CSl*6uMl IM 
CST3*CSl*SUMlKL 
cst4*sni*s>umUm 
IF  (U1.LT.0.)  GO  TO  200 
Xl  1 1W_«XKA¥o*KST  1-CST21 
XI  ll  1M«XKAY0»*CJ»T3^CSTA1 
GO  TO  300 

200  x 1 1 1HL-XKAY0* « CbT  1 *CST2  1 42.*  « CS 1- 1 .4XKYOSO*SUM) 

Xl  1 1 IM»XKAY0#*CaT3-CST*>-SNl-SNl 

300  CONTINUE. 

CST1«<X0/H>+1« 

CST2*CSTI-2« 

xxi irl*csti*csi-xi i me 

XKl  1 lM*-CSTl*iNl-Xl  1 1 1M 
XK12RL*CST2*CS2*XI 12RL 
XKi21M«-CST2*SN2+Xl 12 IM 
bUMIRL'XKl IRL+XX12RC 
SUM  1 1M>XK1 1IM+XX121M 
XK»EAC-.S*<Cb»SUMlHC+SN«SUMl  IM) 

XK 1 MAC* .5*  ( CS»SUM  1 1 M- SN»SUM  l RL  » 

RETURN 

STiAUY  FORM  UF  KtKNtL  FUNCTION 

400  XKBEAL*XO/R 
XK1MAC*0* 

RETURN 

SPECIAL  FORM  OF  XtUNtC  FUNCTION  AT  Y0*0 

SOO  XXREAL-CS 
XK1MAG*-SN 
RETURN 
ENO 


b 059 
6 060 
b 061 
b 062 
d 063 
b 06A 
8 063 
B 066 
B 067 
8 06B 
B 069 
o 070 
b 071 
B 072 
B 073 
d 074 
B 075 
d 076 
b 077 
d 07B 
8 079 
B 080 
B Odl 
b 0b2 
8 0B3 
B Od4 
d 065 
8 0B6 
B 087 
8 086 
B 069 
8 090 
B 091 
B 092 
b 093 
o 094 


SUBNOUTlNE  LGSPANl  X « Y *X£L * YEL  »E  l NT  • W *MAXlNT  *►■  > 


C 

c 

c 

c 

c 


INTEGRATION  OVER  HALF  Of  THE  SlNGULAH  STRIP 
* 1 TH  CAUCHY5  PRINCIPAL  VALUE 

USING  A SIXTH  UEGRtt.  lAGRANGIAN  INTERPOLATION 


NEAL  N 

COMPLEX  XX6AR»F*SUM 
DIMENSION  H4l«iUM(4»«M4) 

DIMENSION  COEF  I A)  «S<4»  *XEL<4  > * Y£L« 4 > *E I NT ( MAX iNT > **< 


DATA  coeF/* i3**72*4*95*-6*a/ 

DATA  S/1 • • *333333 33 •- *33333333* - 1 */ 
SFCNl  A.tt)  ■•25*1  1«4-Ai*<  l •+ai 


MAX 


1NT  ) 


C TtST  TO  SEE  IF  S 1 NGUcAR l T Y IS  TO  LEFT  ON  NIGHT 
C 

IFl1P*0 

IF  ( ( Y4Y-YEL III- YEL ( 2 ) ) • GT • 0 ) IFLIP*I 

EPS»ABSCYEL<2)-YELI 1 ) ) 

DO  10  11*1*4 
10  FI l I »*CMPLX(0**0*» 


c 

c 

c 


ueCIN  SFANWISE  INTEGRATION 
OO  1000  1*1*4 


c wGMPLTE  the  upper  ANO  LO*ER  LIMITS  FOR  THE  CHONDAUSe  INTEGRATION 
C NeVENSE  SIGN  IF  SINGULARITY  IS  TO  THE  NIGHT 
C 

21-SI 1 ) 

IF  (IFl1P*EO*1)  U—Zl 
CST 1 * *5*1 l*-21) 

CST2**5*(  1 • ♦21 I 
A*CSTl*XEL<4l+CSTE*XEL«  3) 

B*CST 1 AXEL 1 l | +CST  £*XEL 1 2 ) 

C|a*S*(U'A> 

C2a*S*(B4A) 

YO* Y-CST l * YEL II I -CST2AYEL 1 2 1 
00  20  11*1*4 
20  SUM! I I )*CMPUX10* *0*  ) 

C 

C BEGIN  CHORD* l SE  INTEGRATION 
C 

00  100  L* 1 .MAX  I NT 
22*E INT (L ) 

X0*X- <Cl*22*C2) 

CALL  XERNEL ( XU  » Y 0 » XkNEAL  * XX 1 MAG  > 

XKBAN  *CMPL  X ( X AR  C AL  * XX 1 MAG  > 

N«  1 I * SFCN ( -2 1 *-4.2) 

N<  2 1 »SFCN( 2 1 *-221 
N(  3 > * SFCNl 2i *221 
N«4)*SFCN«-2i»221 
00  30  11*1*4 

30  SUM!  1 I >*SUM<  I 1 )A*ML)*NC  1 1 )*XXbAR 
100  CONTINUE 
C 

C ENO  CHOHO»ISe  INTEGRATION 

c 


c 

001 

c 

002 

c 

003 

c 

004 

c 

005 

c 

006 

c 

007 

c 

006 

c 

009 

c 

010 

c 

01  1 

c 

012 

c 

S 1 3 

c 

014 

c 

015 

C 

016 

c 

017 

c 

016 

c 

019 

C 

C20 

c 

021 

C 

022 

c 

023 

c 

024 

C 025 
C 026 
C 027 

c 02b 

C 02V 
C 030 
C 031 
C 032 
C 033 
C 034 
C 035 
C 036 
C 037 
C 03b 
C 03V 
C 040 
C 041 
C 042 
C 043 
C 044 
C 04b 
C 046 
C 047 
C 04b 
C 04V 
C C50 
C 051 
C 052 
C 0a3 
C 034 
C OSS 
C 036 
C 0»7 
C 065 


non 


OST«Cl*COeF<  II  c OS® 

Ow  40  11*1*4  C 360 

40  r< i i»«f< n»4.csr*iUM( i n c 06i 

1000  CONTI  Nut  c 062 

00  90  11*1*4  c 063 

90  F(  | I liFIII  1/tPS  c 064 

C 069 

fcNO  SPANWlse  INTEGRATION  c 066 

C 06? 

Uc-TURN  c 066 

END  c 06W 
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I 


\! 

p; 

5 


i jjj 


■ 

% 


bOBHOUT  1NE  P0u  1 uii ( X « Y « XEl. « YEL • fc  1NT  » ft  *MAX  1 NT • )■  • l ANOb-c. ) 

C 

C INTEGHaTION  OVER  A KE&ULAW  TW1ANGULAW  ON  UUAjw I Lt* TiRAL.  tLtMtNT 
C AS  INDICATED  6Y  THE  VALUE  OF  I ANCLE 
C 

NEAL  N 

COMPLEX  X.<BaK  «F  ♦ SUM 

DIMENSION  F( A) »SUM(4) «N( A ) *XEL I A ) • YEL  (4  > 

DIMENSION  ElNT(rtAXlNT) *6(MAX1NT) 

SFCMA.6)».£3*(  1.4A)*( l*4b) 

IFlAC-I ancle-2 

00  100  ll»l«4 
100  Ft  1 1 ) aCMPLX ( 0*  »0*  ) 

X3l*XEL(3)-XELU) 

Y3l * YEL 1 3 )— YEL ( i 1 
CO  TO  (200*300) * 1FL AG 
C 

c triangle 
c 

200  X23»XEL(2)-XEL( 3) 

Y23»YEL(2)-YEU  3) 

c 

C • ACOdlAN  / 4 
r 

XJCbN- .25* ( Y3 1*A23-X3 1 *Y23 ) 

00  2 SO  J»1«mAXIMT 
XJaElNT(J) 

N(2)-.S*( l*+XJ) 

C I a*S*( I •— XJ ) 

WJaCI*W( J> 

DO  210  11*1*3 
2)0  SUM! 11) *CMPLX ( 0*  *0* ) 

OO  230  I*  l *MAX1NT 
XlaElNTI 1 ) 

Nil) ■ »5*Cl • ( 1 *+X 1 > 

N(3)»C1-N( 1 ) 

2ETAaN( 1 ) *XEL ( 1 ) +N( 2 ) *XEL( 2 ) 4N( 3 )*XEL ( 3 ) 

ETAaNt 1 )*YEL< 1 )+M2  )*Y£L  ( 2 ) +N(  3)*YEL(  3> 

XOaX-ZETA 

VOaY-ETA 

CALL  KERNEL cX0»Y0*XKREAL*XK1 MAG) 

XKBARaCMPLX(XXR£AL«XK 1MAG) 

MlaWt 1 )/( YO*Y0) 

00  220  11*1*3 

220  SUM(ll)aSUMim+WI*N(  U**XKBAH 
230  CONTINUE 

00  240  11*1*3 

240  Fill) »F (11)4  YJ*SUM( 11) 

2s0  continue 

GO  TO  400 

c 

c uUAOW ilatehal 
c 

300  X42BXEL (4 ) — XEL ( 2 ) 

Y42» YEL ( 4 ) “YElI 2 ) 

C 

C JACOB 1 AN  / B 
C 


J 00) 

V 002 

J 003 
3 004 
J 00s. 

J 006 
3 00? 

3 OOP 
3 000 
3 0)0 
J C) ) 

3 0)2 
3 0)3 
3 0)4 
u 013 
3 0)6 

3 0)? 

O 0)6 
3 010 
3 020 
3 021 
3 02c 
3 C2J 
3 024 
O 025 
0 026 
u 027 
O 026 
J 029 
0 03w 
3 03) 

3 032 
3 033 
O 034 
0 035 
3 036 
O 037 
0 036 
3 039 
J 040 
3 041 
3 042 
3 043 
0 044 
O 045 
3 046 
3 047 
3 046 
3 04V 
3 OaO 
3 05) 
0 052 
3 053 
O 054 
0 055 
O 066 
0 0S7 
U 056 


J 


f 


I 


1 


✓ 
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XJCbN-.125*lX31*Y42-Y3l*X42» 

l)U  350  I«1 .MAXINI 


DO  310  K«1 *4 
310  aUMlK»*CMPLXl0«»0«> 
ZlaEINTi 1 > 

DO  330  JaltMAXINT 


Z2aElNT«J> 

NC i )>SFCN(-Zl*-Z2) 

N«2»*SFCN1Z1 »-Z21 
N<3)>SFCN(Z1 *12 > 

N(4l>SFCNI-Zl*Z2l  , ... 

^tTA«M  1 »*XELl  1 )4'N*2>*Xfc4-12)4NI  31*XfcL.(3»4Nl4  )4XtH4  » 
ETA«N(  1 »*VfeL«  1 J4Ml2»*VEH2)4'Nl3)*Yfel.«3»4NU»4YEU<4» 


XOaX-ZETA 

YOaY-ETA 

CALL  KERNEL.  < XO  • YO  • XKHEAl. . XK 1 MAG  ) 
XKBARaCMPLX ( XkRlAl • XX 1MAG > 
tfjatft  J)/» YO*Y0) 

DO  320  K«1 «4 

320  SUM«Kl»SUMlKl-M(JAN«K)4XKBAH 
330  CONTINUE 

DO  340  *«1 *4 

340  F«Kl«F(Kl+»(l|AiUMlK> 

3b0  CONTINUE 
400  CONTINUE 

DO  500  *>1 »4 
500  F 1 K 1 «X JCBN*F 1 A I 
RETURN 


end 


J 059 
J 060 
O 061 
D 062 
D 063 
D 064 
J 065 
D 066 
D 067 
O 068 
0 069 
D 070 
D 071 
D 072 
O 073 
O 074 
D 075 
D 076 
0 077 
D 078 
D 079 
D 080 
D 081 
D 082 
D 063 
D 084 
D 085 
J 086 
0 067 


SUBROUTINE  OlU(X«Y«  INFO*XLlFT*FF.UH*AL*UlMAG.NH.NELfcM.NMOOfc*MX» 

C 

C THIS  SuBOUTlNE  COMPUTES  THE  GENERALIZE  FORCt  LOEhF IClcNTS 
C FOR  A GIVEN  LOAo  UlSTRIBUTlON 
C 

REAL  N 

COMPLEX  XL1FT 

DIMENSION  X(NP»  *Y(NP)  « 1NF0(N£LEM*MX) * XL  I FT (NP  tNMODE ) * FF I NP • NMODE > 
DIMENSION  XEL* 41  *YfcL  * 41  *Ft4*7>  • SUM<4«7>  *FI  < 7 ) *N<4  ) *E  INT  <2  ) • « ( 2 I • 

1 UREAL4 7*7) .Q1MAGC7.7) 

COMMON  XM*XK  * BE T AaO tbETA 
DATA  MAXI NT/2/ 

DATA  El  NT/  • S7  7 3 50269  1 UV6ct>  • - • 577  350269  1 t)Vo  tb/ 

DATA  W/2*l ■/ 

SFCN(A«b)a«2S*( 1«4A)**1«4B» 


C 

c 

c 


INITIALIZE 

DO  2000  I-I.NMOOE 
DO  2200  J* 1 « NMODE 
OREALI ItJlBO* 

2200  vilMAGI  1 *J>»0* 

DO  2400  K* 1 « NP 
2400  FF|K.I|«0* 

2000  CONTINUE 

loop  through  the  llements 

DO  4000  L"1 *NLLEM 
IAnGLE»4 

IF  I |NF0(L*4> *E0*0>  I ANGLE >3 
DO  3000  KM.IANGLE 
Kla INFO(L*K) 

XEL (XI »XI  K 1 ) 

YtL<F)»V(XI J 

3000  CONTINUE 

iflag«i angle-2 


c 

c 

c 


% 


E 001 
E 002 
E 003 

004 

005 

006 
007 
006 

009 

010 

011 

012 

013 

014 

015 

016 

017 

018 

019 

020 
021 
022 

023 

024 

025 

026 


E 

E 

E 

E 

E 

E 

E 

t 

E 

E 

E 

E 

E 

£ 

E 

E 

E 

E 

E 

E 

E 

E 

E 


% 


C 

C 

c 


c 

c 

c 


c 

c 

c 


CLEAR  F BEFORE  ACCUMULATION 

DO  100  I J»  l «4 
DO  100  M00E»1*NM00E 
100  F ( | 1 *MCOE>»0* 

X3l »XEL (3  > — XEL ( 1 1 
Y3l»YEL*3)-VEL* 1 I 

BEGIN  INTEGRATION  — OUTER  LOOP 

OO  290  U»l • MAX I NT 
GO  TO  (251*2521* IFLAG 

TRIANGULAR  ELEMENT 

251  XJaEINT ( J) 

N(2) «.5*( 1 • ♦XJ| 

C1-.5** 1.-XJ) 

WJbC1*W( J) 

GO  TO  253 


E 027 
E 02B 
E 029 
E 030 

031 

032 

033 
E 034 
E 035 
E 036 
E 037 
E 038 
E 039 
E 040 
E 041 
E 042 
E 043 
E 044 
E 045 
E 046 
E 047 
E 046 
E 049 
E 050 
E 091 
E 052 
E 093 
E 054 
E 055 
E 056 
E 057 
E 05B 
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no  o 


C UUADR1L.ATESAL  EL  t MEM 

c 

252  22* EINT(J) 

WJ=V»(  J) 

25  3 00  210  11*1*'* 

DO  210  MODE.*  1 * NMODE 
210  SUM  1 1 I .MODE ) *0 • 


C oEGIN  INTEGRATION  — INNER  LOOP 

00  230  I«1  •mAAI.mT 
GO  TO  (231 *23/) • I^LAG 


C 

C 

c 


TH 1 ANGULAR  EuEMEN T 


231 


XI«E1NT( 1 > 

N ( i ) * »5*C l*(l»*Xl> 


N( 3 I =C 1 -N ( 1 I 
XP.N( 1 I*XEL( 
YP»N( 1 I *YEL ( 
GO  TO  233 


1 ) ♦ IN  ( 2>*XEL(2  )+N(3)*Xl  L ( 3) 
1 ) *N ( 2 ) * YEL ( 2 » *N ( 3 > * YEL ( 3 ) 


C 

c 

c 


quadrilateral  element 


232 


2 1 »E 1NT ( 1 > 

N(  1 )»SFCN(-Zl *-i 2» 
N(2>»SFCN(Zl *-Z2> 

Nl 3 ) *SFCN ( Z I .Z2> 

N(4>.SFCN(-Zl«Z2) 

XP»N( 1 1*XEL( 1I+N(2I*XEL« 
YP-N ( 1 1*YEL( 1 )+N(2»*YEL( 


J+N(3)*XtL(3»*N(A)*XtLlA> 

) +N ( 3 > * YEL  <31 +N (A  > *YeL ( A > 


lompute  the  modal  Deflection 


233  FI ( l >■! • 

FI (21-XP 
Fl(3»«XP*XP 
F I ( 4 1 • YP*YP 
F I ( 5 > «F 1 ( 3 > *F  I ( 4 > 
FKNMODE-l  l.YP 

FI ( NMODE >«XP»YP 

DO  220  I l ■ l • I ANGLE 
CST*W( 1 »*N( 1 I ) 

DO  220  MODE. 1. NMODE 
220  SUM(  1 I •MODE)  ■SlMU  I 
230  CONTINUE 


.MOOE>TCST*F 1 (MODE ) 


C 

C 

C 


C.NO  inner  loop 


DO  240  I1«1,IAN0LE 
DO  240  MODE.l.NMOOE 
240  F( 1 1. MODE >.F ( II. MODE) >WJ*SUM( 
250  CONT 1 NUE 


1 1 .MODE) 


C 

C END  OUTER  LOOP 
C 

GO  TO  ( 200 .300 ) . IFLAG 


c.  05V 
E ObQ 
t Obi 


E 062 
E 063 
E 064 
E 065 
E 066 
E 067 
E 066 
t 069 
E 070 
E 071 
E 072 
E 073 
E 074 
E 075 
E 076 
E 077 
E 07b 
E 079 
E OttO 
E Obi 
E 062 
E 0o3 
t.  064 
E Obb 
E 066 
E 067 
E 066 
E 089 
E 090 
E 091 
E 092 
E 093 
E 094 
E 095 
E 096 
E 097 
E 09b 
E 099 
E 100 
E 101 
E 102 
E 103 
E 104 
E 105 
E 106 
E 107 
E 1 Oo 
E 109 
E 1 10 
E 111 
E 1 12 
E 113 
E l 14 
E 115 
E l 16 
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n n n 


T n | ANGt  JL  AW  ELEMENT 

dOO  X23=XEL ( 2 1 "XtL l 3 ) 

¥23  = YtL<2>-YtU  3) 

AH£A..25*< Y3l»X23“X3l*Y23> 

GO  TO  400 

UoAORlLATERAL  ELEMENT 

3G0  X42=XEL< 4 ) -XEL l 2 ) 

Y42=YEL<4)-yEL<2> 

AR£A=. 125*<X3 1*Y42“Y31*X42 

GLOBAL Ut  THE  ELEMENTAL  INTEGRATED  COEFFICIENT  VlLTOW  TO 

t-UHM  THt  INTlORATcD  COEFl-  1C1EnT  vfcX 

400  DO  500  J» 1 • I \NGlE 

K=  INFOIL. J> 

DO  500  MODE«l*NMOOE  _c , 

- , w,  jf\r\c  \ > adr  ( J i MODE  ) 

500  FF ( K.MOOE ) = FF l K .MODE  +AKt« 

4000  CONTINUE 

COMPUTE  THE  GENERALIZED  FORCE  COEFFICIENTS  U«l.0» 

DO  5000  I»1«NMUDE 
DO  SOuO  J*l»NM00t 

DO  5000  K*I,NP  eelx*  I Al ( XL l F T 1 K • J > ) 

UWeALU.J)«Q'JtAL  *»*j  ~ F )#mmag(XL1FT1n.U>> 

dimagi  I «3»  *OlMAU  1 1 »0)  Frnv«»* 

5000  CONTINUE 

C O.V.DE  IMAGINARY  PART  OF  U SV 
c SKIP  the  conversion  for  zero  frequency 

C ,F  IXK.LT. I.E-5)jiC  TO  99V9 

DO  6000  1*1. NMOOe 
DO  6000  J.l.NrtODc 
6000  U l MAG I l * J 1 *Q I MAG  I I »J)/XK 
9999  RETURN 

end 


l 117 
E 1 lb 
L 1 19 
L 120 
L 121 
E 122 
E 123 
E 124 
E 125 
E 126 
E 127 
E 12B 
t 129 
E 130 
E 131 
E 132 
E 133 
E 134 
E 135 
E 136 
E 137 
E 13b 
E 139 
E 140 
c.  141 
E 142 
E 143 
E 144 
E 145 
E 146 
E 147 
E 1 4B 
E 149 
E 150 
E 151 
E 152 
E 153 
£ 154 

E 155 
E 156 
E 157 
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n n n non 


jsUBROuTINl  SCiKMbS  I XNOUE  • YNOUt  i XtL  i Yt.L  «i  l NT  t » t MAX  l NT  i F • 3tL)  ' 001 

c t 002 

C THlb  buBROUT  1 NE  PtRFORMS  THE  INTEGRATION  OF  A RHOMulC  tLEMLNT  ' 003 

C WITH  The  SINGULAR  STRIP  PASSING  THROUGH  lITi-iLR  Tut  LEFT.  » 004 

C MIDDLE  OR  RIGHT  NODe  ' 005 

c * 006 

COMPLEX  F.Fl  r 007 

DIMENSION  XEL<4) «YEL< A) «XEL1 <4 ) « VEL I <4> «F«4 ) .F I (4 ) 1 00B 

DIMENSION  E1NKMAXINTI«WIMAXINT)  ' OOW 

010 

INTEGRATE  THt  RIGHT  HALF  TRIANGULAR  REulON  Oil 

* 012 

CALL  SI  NGUL  ( XNOOfc  * VNODE  » XEL  * YEk.  «E  INT  » W .MAX  I NT  *F  .DEL  ) > OIJ 

f 014 

INTEGRATE  THc.  LEFT  HALF  TRIANGULAR  RtGlON  f Oita 

* 016 

XEl 1(1) SXEL ( 1 ) * 017 

VEL I < I >*YEL< 1)  * 01U 

XEl 1 < 2 ) = XEL ( 3 ) f 0 19 

YEL  I < 2 ) ■ YEL  ( 3 ) ► 020 

XEL 1(3) “XEL ( 4 ) f 021 

YElI <3>“VEL<4»  f 022 

call  Si NGUL ( XHODt » YNODE.  XtL 1 » YLL I «E INT  »W .MAX INT  »F 1 • DEL  I f 023 

pi i )*f< i >+fi m F 024 

F l 3 1 *F  ( 3 ) ♦F  | ( 2 ) F 02E> 

F ( 4 I »F 1(3)  f 026 

RETURN  F 027 

ENO  f 02U 
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SUBROUTINE  SGTRGLtXNODE.YNODE.XEL.  YEL.21NT  . W.MAX 1NT .F .DEL > 

C INTEGRATION  OVER  A utNtKAL  TRIANGULAR  C.LEMENT  WITH  THE 
C SINGULAR  ST«1P  PASSING  THROUGH  EITHER  ONc.  Oh  THu  VERTICES 

C 

COMPLEX  F.Fl 

DIMENSION  XEL l 4 ) » YEL t 4 ) • XEL1 <4  J « YEL lt4>»FC4>«Fl<4> 
DIMENSION  E1NTIMAX1NT) .WtMAXlNT) 

DEL  1 « t YEL t 1 ) - YELl 3 > ) / ( YEL  t 2 ) “YEL  t 3 > 1 
IF  (DEL1.LT. 1.E-5I  GO  TO  1000 


ON » 1 .-DELI 
IF  tDM.LT.l.E- 


5)  GO  TO  1000 


C 

c 

c 


c 

c 

c 


INTEGRATE  THE  LEFT  SIDE  OF  THE  TRIANGLE 

XELl t 1 »«XEL( 1 ) 

YEL 1 < 1 I K YEL ( 1 I 

XELl <2)*0EL1 *XEL(2>4DM*XEL«3) 

YELl (21«DEL1*YEL(2)4DM*YEL<3> 

XELl <3»aXEL<3> 

YELl <31*YEL( 3) 

CALL  SINGUv.  ( Xi.OOt « YNODE  . XEL  1 . YEL  1 »E  INT  • 4 .MAX  1 NT  .F  1 «DEl  ) 
Ft 1 1 »F1 < 1 ) 

F<21 »DELl»Fi (2) 

F<3>=F1<3I+0M*F1(2» 

integrate  the  RIGHT  SIDE  of  the  triangle 

XEL  1 t 3)*XEH  ( 2 I 
YELl  13I*YEL1 I 21 
XELl  < 2 1 “XELt 2 I 
YELl  <2>«YEH2I 

CALL  S INGUL ( XHJDE » YNODE . XEl 1 <YEl1 «e 1 NT  . W .MAX  I NT .F 1»DEL> 
Ft  1 )«FI 1)+Fl l U 
F < 2 1 *F < 2 > +F 1 I 2 I +DEL |*Fl < 3 ) 

Ft  3)aF<3)4DM*Fl( 31 
RETURN 


C 

C DEGENERATE  case  WHERE  THE  TRIANGULAR  element  has  ONE  side 
c COINCIDE  with  The  singular  strip 

1000  CALL  oinguhxnoge«ynode*xel«yel.eint.w.maxinT.f.del> 
RETURN 

end 


COPY  AVAILABLE  TO  DOG  DOES  NOT 
PERMIT  FOLLY  LEGIBLE  PRODUCTION 


G 001 
G 002 
G 003 
G 004 
G 00b 
G 006 
G 007 
G 006 
G 009 
G 010 
G Oil 
G 012 
G 013 
G 014 
G 01b 
G 016 
G 017 
G 016 
G 019 
G 020 
G 021 
G 022 
G 023 
G 024 
G 025 
G 026 
G 027 
G 026 
G 029 
G 030 
G 031 
G 032 
G 033 
G 034 
G 035 
G 036 
G 037 
G 036 
G 039 
G 040 
G 041 
G 042 
G 043 
G 044 
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non 


bUQWOUT  IN-.  Si  NGUL  ( XNODE . YNUDE  . XEL  . YEL  . E 1NT  . W .MAX  1 NT  . F . JtL) 


C 

c 

c 

c 

c 


INTEGRATION  over  the  special  type  of  triangular  element 
tlY  A VERTICAL  LlNt  AND  TWO  OTHER  STRAIGHT  LINES 
THE  TRIANGLE  may  BE  SINGULAR  or  REGULAR 


lounged 


complex  F • F 1 
DIMENSION  F(4I  |FU4I 
DIMENSION  XELI  A)  tYEL<  41  «XEL1  14  > 
HI A.U)*OcL*A+CM*B 


YEL 114) »t INTI  MAXI  NT ) • W ( MAX  1 NT  1 


CM=1.-DEL 

IF  I ABS( YEH 1 )-YEL(2>  l.GT.l.E-51  GO  TO  10 
IF  ( ABSl YNOOE-YEl l 3) 1 «LT» 1 *E-b 1 GO  TO  1 
IF  ( AOS<  YNODE~YtH2»  1 «LT»  1 «E-b  > GO  TO  3 
10  CONTINUE 

IF  1 ABSl YNODE- Y tL  12  1 1 »LT» 1 «E*S • GO  TO  2 
IF  ( ABSl  YNOdE“TEL«3)  1 <LTt  1 *£“b)  GO  TO  4 


NON-SInGULAR  TRIAimGlE 


call  POLYGNl  XilODt  .YNODE  • XEL*  Yel.  «E 1NT  « W « MAXI  NT  «F  .3) 
RETURN 


C 

C 

C 

C 


triangle  pointing  left  with  singular  strip  passing  through 
THE  VERTICAL  SIDE 


1 


XtL' < 1 > *XEL ( 3 ) 

YElI‘  1 )=YEH3) 

XElI 121  =HI XEC I 11  »XEL«3>  » 

YliL  K2  >*Hl  YE’_  ( 1 ) « YEL(  3>  > 

XELI  «3»*H«XEL<  21  «XEH3>  > 

YEL l l 31 *H l YELl  2 1 t YlLI 31 > 

XEL  l (4I*XEH31 

YELl  (41*YtM31  „ ...to. 

CALL  LG^.PAN  ( XNOOfc  • VNODE  . XEL  1 • *tL  1 *E  1NT  . w .MAXl  NT  .F  l 1 

FI 1 1 *OEL*F 1(2) 

F 1 2 • *DEL*F l (3) 


F l 3 > =F  1 l 1 > 4-F 1 C 4 1 +CM#  (FH21+FK31) 

IF  1 ABSl 1 . -DEL) «LT. 1 »E-5)  RETURN 
XtL 111) =XELl ( 2 1 
YELl l l Is YELl ( 2 1 
XEL 1141 =XEL l ( 3 ) 

YELU41-YEL1  ( 3) 

XELI  121*XEH  1) 

YLL  l 12  l=YEH  1) 

XEL  l 131 =XEL ( 2 1 

CAll  POL YGN ( XGDOt . YNODE . XEL 1 . YtL 1 « E INT • W .MAX 1 NT . F 1 . 4 1 
F 1 1 1 *F 1 1 1 fF 1 C 2 > +DEL*F 111) 

F12) *F121+Fl 1 3 1 ♦ JEL*F 114) 

F l 3 l =F 1 3 1 +CM* (Fill) +F 1(411 


RETURN 


c TRIANGLE  pointing  wight  with  singular  strip  parsing  through 
c THE  VERTICAL  Slot 

c 

2 XEl 1 l l 1 *Hl XELI 1 1 .XELI 2 1 > 

YEl  1 l 1 >=Hl  YEU  1 1 »YELl  21  t 


H 001 
H 002 
H 003 
H 004 
H 00b 
H 006 
H 007 
H 006 
H 009 
H 010 
M C 1 1 
H 012 
H 013 
H 014 
n 01b 
H 0 16 
M 0 17 
H 016 
H 019 
H 020 
H 021 
H 022 
h 023 
M 024 
H 025 
H 026 
H 027 
h 026 
H 029 
H 030 
H 031 
H 032 
h 033 
h 034 
R 035 
n 036 
H 037 
H 038 
H 039 
H 040 
H 041 
H 042 
H 043 
H 044 
H 04b 
H 046 
H 047 
H 046 
H 049 
H 050 
H 05 1 
H 052 
H 053 
H 054 
H Oob 
h 0b6 
H 057 
H 056 
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COPY  AVAILABLE  TO  DOC  DOES  HOT 
PERMIT  FULLY  LEGIBLE  PRODUCTION 


permit  fully  legible  producth 

XLLI(3l«XtU  21 
YEL 1 < 3 1 ■ YEL ( 2 ) 

XEl»  XEU  31  »XEL<21  1 

VEL  1(41  *H  ( VEL  ( 3 ) * VEL  (2)1  MAV1K,TP,, 

CALL  LG5PAN1  Xi’tODfc.  • YNOJE  • XEL  1 • YtL  1 • 1 1NT  • V »MAX  INT  *F  1 ) 

F l 1>»DEL*F1UI 

K<2»«FJ  (2MFl(3MCM*<Fl(l  UFlHI  » 

F ( 3 1 «DEL*F 1 (4) 

IF  (ABS(  1«-0EU«LT»1»E-51  RETURN 
XEL1(2>*XEli(1) 

YEl 1(21*  VEL 1(1) 

XtLl ( 3»»XELl (4) 

VcLll3l*VELll4) 

xel  i « i i »xen  1 1 

YEl 1 ( 1 >»YEL( 1 I 
XEL  1(41 >XEL ( 3 ) 

call  * POL  YGN«  XIWUL  • VNOOE  « XdL  1 • *tL  t »L  t NT  • * •NaX  INT  .F  l . 4 l 
F< 1 1 «F( It ♦Fl ( 11+JEL4F1 <21 
F (2 ) »F ( 2 ) *CM*( F1I2I+K 1 < 3 » » 

F<3>«F(3> *F i ( 4 1 ♦OEL*F 1(31 
RETURN 

C triangle  pointing  left  with  singular  strip  passing  through 
c the  vertex  OPPOSITE  To  the  vertical  side 
c 

3 XEL 1<1I *H<  XEL I 3 1 fXELI III 
VtLl!n«H»VEL<3>  •VELUM 

X£LII2I«XEH1) 

YEL 1(21 »VEL ( 1 1 
XEl 1(31 “X6L( 2 I 
YEL 1*31 *VEL( 2 1 
XELl  I4MHIXEU  31  *XEL<21  1 

YELH4)»H(YEL(3I  •YELI2I  1 

CAcu  LG5PAN1  Xi'IJOfc  • YNOJE*  XLL 1 « Ycu  MEINT»W»*  4X1NT«F  I I 

F(1)«FK2I*CM*F1UI 

F(21«FI»3MCM*F1  141 

F ( 3 1 »DEL  •(Fl(ll+FH4l  1 

IF  ( ABS<  M-OEE)  •LT«l»E-5»  RETURN 

XELl  (3>*XEL| I 1 1 

YEL  1(3  J * YEH  ( X ) 

XtL I 1 I l*XEL|(4l 
YELI < I »«VEci (4) 

XEL 1(21 *XEL ( 3 i 
YEL  1 ( 2 ) * YEL  ( 3) 

call  POLVGN  (X  NODE  *VI400E*  XEL  I • VELI  *C  INT  » * tNAXlNT  »F  M3! 

F ( M *F  < II ♦CM*F 1 ( 3 1 
F ( 2 1 »F  <21 >CM*F 1 I 11 
F ( 3 1 »F  (31 +F 1 ( 2 MOEL*(Fl  (|HF|(3)I 
RETURN 

C TRIANGLE  POINTING  iUGHT  with  SINGULAR  STRIP  PASSING  THROUGH 
C THE  VERTEX  OPPOSITE  TO  THE  VERTICAL  SIDE 
C 

4 XELl ( 1 1 *XEL ( 1 1 
YtLl ( I l“YEL(l» 


h 05V 
H 000 
H 061 
H 062 
h 063 
H 064 
ri  065 
H 066 
H 067 
H 066 
H 06V 
H 070 
H 071 
H 072 
H 073 
H 074 
H 075 
H 076 
H 077 
H 076 
H 07V 
h 060 
H 061 
H 062 
H 063 
H 064 
H 065 
H 066 
H 067 
H 066 
H 069 
H OVO 
H 091 
H 0V2 
H OVO 
H 094 
rt  0V5 
H 096 
H 097 
H 098 
H 099 
H IOC 
H 101 
H 102 

H 103 
rt  104 
H 105 
H 106 
H 107 
H 108 
H 1 09 
H 110 
H III 
H 112 
H 113 
rt  114 
H 115 

H 116 
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- 


H 

117 

XtLl  12»=H<XEU2>  *XtL«  1 ) > 

H 

1 ia 

VtLl  «2»=H«  YeH  2>  »YE.Li  1 > ) 

M 

119 

XtLl ( 3) »H«  XEL<  2 > »XtL»  3) 1 

H 

120 

YtLl  »31»H1  YEH  2)  tYELl  3)  > 

M 

121 

XELl (4)*X£L| 31 

H 

122 

CALL*  LtoPANlXlloUti*  VNODE.  X£L1.  VELl.t  1NT.  W.MAX1NT.F  1 > 

H 

M 

123 

124 

F«1  )*F11 1 >*CM*Fii21 

H 

1 25 

Fig) »DtL*  1 F 1 ( 2 1 +F  1 1 3 1 1 

H 

120 

F« 3) »F 1 1 A »+CM*F 1( 3) 

H 

127 

IF  ( AbSU  .-Dtl- ) «i-T.  1 «E-5>  NETUWN 

H 

1 2t> 

XELl 1 1 )=XLL1 (2) 

H 

129 

VELl«l)=VEt-l<2) 

H 

130 

XtLl <2 1 *XEH  2 • 

rl 

131 

CALL  * POtYGN  < ANODE  . VNODE  . XEL  l » VtL  1 » 1 1 NT  . « . MAX  1 NT  .F  1 . 3 > 

M 

M 

132 

133 

F(ll.F(l)+CM«Hlll 

M 

134 

F 1 2 > *F  ( 2 > +F 1 ( 2 )+DEL*<F  H1HF1I3I) 

rl 

135 

F * 3 ) *F ( J ) +CM*F 1 ( J» 

H 

13b 

WE TOWN 

H 

137 

END 

i 


SUbKOUT I Nt  T ADLC  1 COEF .DEL T A. t I N T • * • "*X l NT .DEL . 1 MAX. JMAX .LMAX ) 

TrilS  SUBROUTINE  CKEATcS  A TABLE  Of  THE  WEIGHTED  *t-RNLL 
LUt-FFlCltKTS  FOR  A Kc^ULAP  CHARACTERISTIC  MLSH 

^CHElImENT  lT  uil-U^Y  UEF  l NED  bY  A PAIR  UF  REcAT, VE  INDICES 
rut  L.,(-iu  )F  THt  SIDE  OF  THE  CHARACTERISTIC  eLEMeNT 
mAx  1 5 J T HE  MAXIMUM  POSSIBLE  NUMBER  UF  ELEMeNTS  IN  THe  CHORDjIISe  OJPt-C 
22  »1  the  MAX, MUM  POSSIBLE  NUMBER  OF  ELEMENTS  IN  THt  MACH  LINE  UlRLC 
LMAX  IS  THE  MAXIMUM  PcSSlULt  ENTRIES  TO  The  Table 
The  TABuE  lb  CumPaCTlV  STuREO  INTu  A mECTANcweAR  AmkAY 

COMPLEX  F « COEF 

DIMENSION  F < 4 > • LdE> ( A .LMAX ) 

DIMENSION  XEL 14, « YLL14, it l NT l MAX l NT  » *4 (MAX l NT  ) 

COMMON  XM.XK ♦BfcTAbL.BETA 

LINEAR, I.U»«1 l l-l »*, UMAX* JMAX-l  » »/2  + U 

YOEL*DELTA/XM 

XDEL«bETA*YOEL 

DO  20  LL* l • l MAX 

DO  20  MM«lL,JMAX 

; utMNE  the  element  EUL-TlON  aeCORDINU  To  ITS  «eeaTIVe  INDICES 

l_*L  l near  (EE,  Mi-1) 

XeE<  1 I ■ ( 2“EL“MMl *xdel 
XEL ( 2 i * XEe ( l ,-XUtL 
XEl.  ( 3 ) ■ XfcE  , ? ,-XUEL 
XEL <4, -XEL, £ I 
YEL  < l 1«(MM-lLI»YUEL 
YEL( 21«YEL , l I ♦ YUEL 
YEL(3I»YEL( l ) 

YEL ,41* YEL ( 1 l -YUEL 


CHARACTERISTIC  ELtMtNT  CDT  bY  Thfc  S , NGUlAR  STRIP 


lr  ((MM-EE).LE.l,  LALL  SGRHbS , XNODE  » YNOUE  • AEL  • YEL  «E  INT  • 4 .MAX  1 NT  • 
t F , dEl ) 


HeUDlAR  charac t ek i s T l c element 


IF  , (MM-LL,  .GT«H  CALL  POLYON,  XNODE.  YNODE » *EL • YEL »t INT • 4 .MAX  1 NT • 
( F ,4  , 

DO  20  l * l » 4 
COEF, I »L , *F , I , 


20  CONTINUE 
RETURN 
END 


Smluif^JU110  ms  H0T 

PERMIT  FULLY  LEGIBLE  PRODUCTION 


i 001 
i 002 
I 003 
1 004 
1 005 
1 006 
l 007 
l 00b 
I 009 
1 01U 

l Oil 

; 012 

1 013 
1 014 

' CIS 
016 
. 017 

. 0 1 B 

014 
020 
021 
0 22 
023 
1 024 

I 025 
1 026 
l 027 
1 02B 
l 029 
1 030 
* 031 
> 032 

033 
I 0 j4 
i 035 

036 

037 

03o 

OjV 

040 

041 

042 

043 

044 
. 045 
1 046 
1 047 
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C32=COb< A«ul  5 057 
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SPECIAL  OF  fctkwEl-  FUNCTION  AT  YO*0 


